 function [nORr] = plotEntropyFix( dirName ,fileNum,M,N,dataMode)
% determine the optimal parameter between M and N, by maximal entropy principle
% dataMode = ...
% 1: normal data
% 2: flow field data
%---------------------------Preparation-------------------------------%
global topological
[position,velocity,num]=readText( dirName,fileNum,dataMode);

S = sum(velocity) / num; % now S is the average velocity vector 
SL = norm(S); % SL equals the model of S
unitVec = S/SL; % row vector
sL = velocity * unitVec'; % projection of velocity, column vector
pai = velocity - sL * unitVec ;  % get the longitudinal velocity "sL" and the perpendicular velocity "pai" 

distance = zeros( num,num );
expCorr = zeros(num,num);
for i= 1 : num
    for j = 1: num
        distance(i,j) = norm(position(i,:)-position(j,:));
        expCorr(i,j) = sum(velocity(i,:).*velocity(j,:));
    end
end % acquire distance and correlation between i and j, programme can be improved.

[logicI,logicB] = distinguishBoundary(position,distance);
% logicB = zeros(1,num) == ones(1,num);
% logicB(1:100) = true;
% logicI = ~logicB;
numI = sum(logicI);

%--------------------------Compute Entropy---------------------------%
if topological == 1
    step = 1;
elseif topological == 0
    step = (N-M)/20;
end % "step" depend on situation: neighbor num or visible radius
im = (N-M) / step +1; % test times from M to N

entropyArray = zeros(1,im);
ncAver = zeros(1,im); % export the average of nc for comparing with topological assumption
for i = 1 : im 
    n_r = M + (i-1) * step; % varible "nORr"
    [logic, ~, corrSum1] = computeIntCorrelation( num ,n_r ,distance ,expCorr );
    ncAver(i) = (sum(sum(logic))/num);

    N_OO = 0.5 * (logic + logic'); % when i = j ,logic(i,j) = False
    N_IB = N_OO(logicI,logicB);
    N_II = N_OO(logicI,logicI);
    AdJ_II = diag(sum(N_II,2) + N_IB * sL(logicB)) - N_II;
    %h_I = N_IB * velocity(logicB,:); % column vector
    hL_I = N_IB * sL(logicB);
    hP_I = N_IB * pai(logicB,:);
    J =(numI-1)*( 0.5 *sum(sum( AdJ_II^-1 .*( hP_I * hP_I'))) - 0.5 * (norm(sum(pai(logicB,:))+...
        sum( AdJ_II^-1 * hP_I)))^2 / (sum(sum(AdJ_II^-1))) + sum(hL_I) + ...
        0.5 * sum(sum(N_OO(logicB,logicB).*(velocity(logicB,:) * (velocity(logicB,:))'))) +...
        0.5 * sum(sum(N_II))  - 0.5 * sum(corrSum1 .* sum(logic)) )^(-1);
    % % doubt?
    % disp( 0.5 *sum(sum( AdJ_II^-1 .*( hP_I * hP_I'))) );
    % disp( 0.5 * (norm(sum(pai(logicB,:))+sum( AdJ_II^-1 * hP_I)))^2 / (sum(sum(AdJ_II^-1))) );
    % disp( sum(hL_I) );
    % disp( 0.5 * sum(sum(N_OO(logicB,logicB).*(velocity(logicB,:) * (velocity(logicB,:))'))) );
    % disp( 0.5 * n_r * (numI - num*intCorrelationI) );
    % disp('J=');disp(J);
    % disp(J*n_r);

    A_II = AdJ_II * J;
    A_II_AdjustNum = max(max(A_II))/2;
    A_II_adjust = A_II/A_II_AdjustNum; % adjust to avoid overlarge A_II
    lnZ = 0.5*J^2* sum(sum(A_II^-1 .*(hP_I * hP_I'))) - log(det(A_II_adjust))-numI*log(A_II_AdjustNum)-log(sum(sum(A_II^-1)))-...
        0.5* (norm(sum(pai(logicB,:))+J*sum( A_II^-1 * hP_I)))^2 / (sum(sum(A_II^-1)))+0.5*J*sum(sum(N_II))+...
        J*sum(hL_I) + 0.5 * J*sum(sum(N_OO(logicB,logicB).*(velocity(logicB,:) * (velocity(logicB,:))')));

    % disp( 0.5*J^2* sum(sum(A_II^-1 .*(hP_I * hP_I'))) );
    % disp( log(det(A_II_adjust))+ numI*log(10) );
    % disp( log(sum(sum(A_II^-1))) );
    % disp( 0.5* (norm(sum(pai(logicB,:))+J*sum( A_II^-1 * hP_I)))^2 / (sum(sum(A_II^-1))) );
    % disp( 0.5*J*sum(sum(N_II)) );
    % disp( J*sum(hL_I) );
    % disp( J*sum(sum(N_OO(logicB,logicB).*(velocity(logicB,:) * (velocity(logicB,:))')))  );
    % disp('lnZ=');disp(lnZ);

    entropy = ( -lnZ + 0.5* J * sum(corrSum1 .* sum(logic)) )/num;
    entropyArray(i) = entropy; 
end

nm = find(entropyArray == max(entropyArray));
nORr = M + ( nm-1 )*step;
ncAverage = ncAver( nm );

%--------------------------Plot----------------------------------%

figure(4)
plot(M:step:N,entropyArray)
if topological == 1
    title('entropy with nc')
elseif topological == 0
    title('entropy with visible r')
    text(M+5*step,min(entropyArray),['ncAverage=',num2str(ncAverage)])
end
text(M+step,min(entropyArray),num2str(nORr))
end


